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ABSTRACT 


Tn this investigation minimum weight design of general wing 
structures satisfying strength, stability and frequency requirements 
is attempted. A computer programme package is developed for the 
automated minimum weight design of a ^^rametric wing structure satis- 
fying the strength, stability, frequency ard flutter requirements. 
Attempts are made to reduce the number of design variables, \4iich 
in turn reduce the required computational effort, by the help of a 
parametric study. This study indicates that uniform mass distribution 
can be considered for the spars and ri-bs without ary appreciable 
change in the behaviour constraints. A typical vang structure is 


optimiz'sd satisfying the above requirements^ The feasibility of 
employing linearly-approximated r e dosignB is investigated. Off-design 
charts are obtained by a parsmetric study about the final optimum 


design point. Further improvements that can be incoipjorated in this 


study are also suggested. 
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CHAPTER 1 


INTROnJCTIOH 

During the design stages of an aircraft, it is customary to 
base the structural design primarily upon the requirements of 
strength and stiffness. The dynamic and flutter requirements are 
then satisfied by making modifications to the strength design by a 
process of trial and error. This procedure is adopted mainly 
because of the lack of a suitable design procedure that can handle 
all the requiirements simultaneously and also the design modifications 
to satisfy the flutter requirements are not extensive for subsonic 
aircraft. However, in the design of large, hi^ perfomance aircraft 
dynamic and aero elastic requirements pl^ nearly as prominent a 
role as the strength. Cases can be cited when the required maigins 
on flutter speed coifLd be met only through the addition of several 
thousand pounds of material to a wing which had already met static 
loading design conditions. Also it is a known fact that the 
aerodynamic parameters like sweepback angle, aspect ratio and 
thickness to chord ratio influence -ttie dynamic and aeioelastic 
performance. Hence these parameters can also be included when the 
design modifications are made to meet the dynamic and aeroelastic 
requi3?Qiients, as long as the resulting changes in the aerodynamic 
parameters does not appreciably change the aerodyn am ic performance. 
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Although techMques for the analytical prediction of dynamic and 
aeroelastio properties of a given structure are hi^ly developed, 
the introduction of such features into fonnal structural procedures 
has not been given much consideration untill recently. 

1 « 1 Review of Literature 

mring the last few years, several structural optimization 

investigations have been reported that deal mainly with one 

particularly troublesome behaviour constraint. Turner^ developed 

a piocedure for determining relative proportions of selected 

elaaents of an aircraft structiore to attain a specified flutter 

speed with mi n im u m total mass, Lagrange multipliers were anployed 

to introduce the constraints and the resulting systan of nonlinear 

equations was solved by the Hewton-Eaphson method to determine the 

masses of the elements of the system, Ihe optimization of one 

dimensional structures like beams under aeroelastio constraints 

2 

was reported by Ashley et al. In this work, the authors found 

the solution by reducing the original problem to a set of first 

order differential equations by the application of Euler- Lagrange 

3 

equations and the transversality conditions- Haftfca replaced a 
continuous flutter constraint by the minimum value of the constraint 
over the range of values of the parameters, thiis reducing the 
computational effort. Approximate expressions are derived for the 
minimum value of the parametric flutter constraint. 
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itoz axid Kapoor reported a capability for minimum weight 
optimum design of planar truss-frame strucimires with distributed 
and concentrated mass. Inequality constraints were placed on the 
maximum dynamic displacements and stresses and natural frequencies 
of the structure were excluded from certain bands, A direct 
Optimization method (the method of feasible directions) which 

consists of a design-analysis cycle was used to solve the optimization 

5 

problem, McOart, et al, developed a steepest-descent boundary 
value method for the design of structures with constraints on strength 
and natural frequencies, A computational algorithm was developed 

6 

which in^ilemented the steepest descent method, Sippel and Warner 
investigated the optimization problem of simple 2 to 3 elaaent 
structures of sandwich type with frequency constraints. They 
considered the axial vibration, fleimral vibration and the rotational 
inertia effects of tip masses in detail, !Eae minimum weight design 

of structures with a specified natural frequency requirem^t was 

. 7 8 9 

considered by Ihrner , Zar^amee and Eubin . 

10 

Kicher reported a capability for findirg the minimum weight 
design of integrally stiffened cylindrical shells by considering 

gross buckling, panel buckling, and skin and rib buckling as behaviour 

11 

constraints. Zarg h amee used a combination of the Hosen gradient 
proj ection method and steepest-descent, ^ternate steps method of 
Schmit in the minimum weight design of a discretized structure with 
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a general stability constraint. Simites anS Unghbhafcorn , 

by judicious choice of the objective ftinction and proper grouping 

of the design parameters accomplished the minimom weight design 

of stiffened cylinders under axial compression in two phases. 

Hie first optimization phases yields design chart Sj which are 
then used in the design phase to arrive at the Tnin-Tmurn 

weight. 

Recently there have been some attempts to include bo1±i 

strength and aero elastic requiroaents in the optimum structural 

13 

design problaas. Schmit and Ihornton considered a highly 
idealized double wedge airfoil with the total propulsive woife as 
the objective for minimization. Ohe root angle of attack, tip 

. V 

deflection, root stresses and bending- torsion flutter Mach number 

were considered as behaviour constraints with airfoil thickness and 

ch^rd length as the two design variables. An automated preliminary 

design procedure for simplified wing structures has been reported 
14 

by Stroud, et al . By idealizing ihe wing as an isopropic sand- 
wich plate with a variable thickness distribution, the authors used 
an interior penally function method to fiM 1die minimum weight 

cover thickness distribution. Strength, stability and flutter 

15 

constraints w^re used, Siles developed a procedure to include 
the interaction of external shape, aerodynamic loads, structural 
geometry and fuel mass in the fully stressed design of thin aero- 
dynsmic surfaces. Eao comidered strength, stability, frequency 
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and flutter constraints under diffe 2 rent flight conditions (varying 
fuel weight). His design parameters include ihe thickness of 
covers, webs and ribs, and flange areas. 

To iEip3X>ve upon the computational effort and expedite convergence 

number of investigators have suggested various refinements to the 

standard algorithms and different foimulations of the basic optimization 

problems. Optimum structural design studies can, with few exceptions, 

be placed into three broad categories. In the first method the stated 

problem is attempted directly by calculus of variations. Although 

this approach is one of formal power and elegence it has been successful 

17 

primarily with plates and portal frames . This method has been 

applied with li m ited success to standard structural systems optimization 

because of the great analytical complexity. Studies in the second 

18 

category are based on Shanley's structural index approach. Ghis 
procedure ha.s b een applied successfully to a variety of minimum weight 
design problems with strength and buckling constraints. This method 
has not been explored in the case of statically indeterminate structures. 
In the third and the most general category optimization of structural 
systems is examined on numerical basis using the tools of linear, 
non-linear and dynamic programming. Ghis approach pioneered by 

19 

Schmit has dealt successfully with a wide range of structural opti- 
mization problems and has been largely responsible for the increasing 
involvement of optimum structural design in traditional design meHiods. 
Although this method is quite successful with a wide range of design 
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problems perticular attention must be given to keep down the compu- 
tational effort and reduce the numerical inaccuracies due to 
instabilii^r while dealing with large structural systems (50 or more 
design variables ), 

. 20 

picket, Jr, et al. developed a procedure to reduce the 

number of design variables in the automated structural ^nthesis 

of large ^sterns. A rational procedure based on the external 

loads and constraints on the system was developed for generating the 

reduced set of coordinates (depend upon the linearly independent 

design like fully stressed design, fully deflected design etc, ). 

Ihus the complete design problem is replaced by much siB5)ler opti- 

21 

mality criteria. Another technique suggested by Pope to reduce the 
computational effort and expedite the convergence, is to cast the 
original mathematical pioblom as a sequence of linear programming 
problems, 

1*2 ScQ-pe of Ihe Present Wock 

As stated above most of the liteiature is ooncerned with the 
optimization of either highly idealized structures, or with structures 
that deal primarily mth one type of h Saviour coiastraint, or 
different formulations of the basic optimization problem to expedite 
the convergence. Tery few investigators considered all the behaviour 
constraints; sifflaltaneously. 
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In the present investigation, the feasibility of including 
multiple behaviour constraints in the optimum design of a symmetric 
wing structure, with design variables including the sweepbadc angle, 
aspect ratio and thickness to chord ratio apart from usual gauge 
parameters is demonstrated, ilso the mass distribution of the wing 
is assumed to decrease linearly towards the tip of the wing, Sie 
inclusion of this feature alone may not accaaplish a fully- stressed 
design, but nevertheless helps in that direction and also some 
saving in the weight can be envisaged over a uniform mass distribution. 
!Ehe gauge parameters, like cover thicknesses and flange areas are 
assumed to vary linearly along the span of the wing. A parametric 
study of the behaviour quantities (maximum deflection, stresses 
induced etc. ) is conducted to study, how they are influenced by the 
design variables, before the actual optimization problem is taken upl 
This study suggested that the n^ss distribution of the span and 
ribs can be taken as uniform throughout the wing, reducing the 
number of gauge parameter variables to 6 in the subsequent optimization 
problaai After the optimum is found, off- design charts are made 
by coaaducting amther parametric study. These charts help to assess 
the penal-ty in the minimm weight whenever the optimum design is 
slightly perturbed. 

A requirement of this nature demands a number of capabilities 
to be assembled, KLrst the structure has to be represented by a 
suitable model to determine the necessary static, dynamic and aeroelastic 
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characteristics. This structural idealization must be accurate 
enough and must be suitable for the embedment into the iterative 
design analysis loop of the optimization procedure. Finally for 
the solution of the mathaaatical programming problem, the partial 
derivatives of the behaviour constraints with respect to the design 
variables have to be detemined. 

In the present work, the finite element method, with 3 
different kinds of elements, is used to model the multiweb wing 
stiucture. The constant stress triangle membrane elements and 
the rectangular shear panels are used to idealize the ^in and web 
respectively. The stringers are represented as axial force 
members. Elastic buckling constraints are introduced by treating 
a typical portion of the wing skin as an isentropic stiffened 
plate. The second order piston theory is used in finding the 
flutter Mach number of the structure. 

The constrained optimization problem is cast as a nonlinear 
mathematical programming problem. The design variables include 
the sweepback angle, aspect ratio, thickness to Chord ratio, 
cross sectional areas of the pin-jointed bar elements, thicknesses of 
spans, ribs, and the two linear constants which represent the skin 
thicknesses, Qhe interior penally function method, with a variable 
metric unconstrained minimization technique, is used in the minimum 
weight design probleam. The minimizing step lengths are determined 
using cubic inteipotation technique, Partial derivatives are 
computed using finite difference foamila. 
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Bae thesis is organized into six chapters. The objective 
function and the design criteria are described in Chapter 2 and the 
problem of minimum wei^t design of wing structure is formulated 
as a nonlinear mathematical programming problaa. 

In Chapter 3 analytical idealization of the wing structure 
is considered. A method of including the elastic buckling constraint 
in the optimization problem is described. [The derivation of the 
necessaiy inertia, stiffness and aerodynamic matrices for the foimu 
lation of flutter equations has also been presented in Chapter 3. 

Chapter 4 deals with the optimization algorithm used in the 
present woifc. Reasons for using the penalty function formulation 
with a variable metric unconstrained minimization method are discussed. 

The results of the present investigation are presented in 
Chapter 5, the conclusions drawn from the present study and the 
recommendations for the future investigation are presented in Chapter 6. 

The derivation of element stiffness and mass matrices for the 
three types of elements used in the present work and the canputer 
programme package developed in the present work are briefly described 
and presented in Appendices A & B respectively. 



CHAPTER 2 


PORMULATION OP 2HE OPTIMIZATION PROBIEM 

Mien a neans of predictiiag the behavioxir of ary design within 
a particular design concept is available, linitations on the perfomance 
and other external constraints on the design can be stated, and an 
acceptance criterion can be established, then it is possible to cast 
the design nodification problea in the fom of a mthanatical progra- 
ming problen, A nathenatical programing problen is one in which 
a nultivariate function f(?) , where X is a n-dinensional design 
vector, is to be optimized subject to given constraints g^(?) £0, 

3 = l,2,««»n. The function f(?) is called the objective function 
and its choice is governed by the nature of the problen. 

2, 1 Objective Punction 

The nature of the structural design problen is such that there 
will usually be nany designs that perform liie specified functional 
purposes adequately. The objective function in structural design 
optimization represents a basis for choice between alternate acceptable 
designs, Wei^t m i n i m ization is often taken as the objective of 
structural design optimization. This is in part due to the fact that 
weight is readily quantified, partly reflects material costs and its 
saving can be converted directly into increased payload or langew 
The minimirium weight of the structure is, therefore, considered as 


the objective function. 
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2 * 2 DesiiSgi Criteria 

Characterization of structural design philosophy imrolves-identi- 
ficatiou of fche kinds of failuro mo das to be guar dad againsto 3?he 
da sxgn cri’t/erxa ijo be satisfied are fiie ^^llowin^ 

(1) She elastoc deflection at tiie tip of the wing must not 
exceed certain prescribed limits^ 

(2) The stresses induced at the root as well as tip of the wing 

should not exceed the yield stress of the material. This constraint 

IS ingjortant since the cover thickness decreases from the root to 
the tip, 

(3) Bhe first few natural frequencies of the stiucture are to 
lie within certain bounds. 

(4) Flutter frequency and Mach number should be away from the 
specified value. However, this constraint is not included in the 
final solution because the double eigenvalue analysis of flutter 

is to be perfoimed n+i times whenever the gradient of the function 
IS computed, which is highly time consuming (IM 7044 computer is 
used). However, the developed pzogramme has the capability of 
including this constraint also, 

2.3 resign Variables 

Hie objective function and the design criteria being known 
lie struotaral optiMzatlea pioblao o* the edag oaa be east as a 
matheaatieal piagreemiag pioblea once the desiga TOidebles are chooser. 
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It is proposed to include sweepback angle, aspect ratio and thick- 
ness to chord ratio as the design variables in order to assess the 
feasibility of optimization in the presence of these variables and 
also to study the inter action of earternal shape and mass distribution 
of the wing, !Ilhe individual thicknesses of the plate elements, 
cross sectional areas of the pin. jointed bars foim the rest of the 
variables. Because of large number of finite elements involved, 
a group of elements are often specified by one design variable. 

One convenient means by which the m design variables of the opti- 
mization problem can be related to the individual thicknesses of 
of plate elements, cross-sectional areas of flanges and weights of 
tuning masses is thiou^ a design correlation table^^. This design 
variable correlation table imposes additional 1 j ritri ng constraints on 
the Optimization problem which can be handled rather easily. Hence 
five pairs of variables, each pair (x. and X. of X.- x X. ) 

1 l+T 1 l+T 

representing the linearly decreasing (in spanwise direction) mass 
distribution of cover plate elements, web elements of spars and ribs, 
and flange elements of spars and ribs respectively, are considered. 

Thus to start with, there are ten variables representing the thick- 
ness of the cover plate elements, web elanents of the spars and 
ribs and flange areas of the spars and ribs. However, these variables 
are further reduced to six, since the preliminaiy parametric study 
showed that the linearly decreasing mass distribution of the spans 
and ribs can be replaced by uniform mass distribution without 
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incurriog aiy appreciable change in the behaviour constraints. 

Thus in the present optimization problCTi nine design variables 
(three variables representing the external shape and six variables 
representing the mass distribution of the wing) are considered. 

2.4 Optimization Problem 


Mathematically the optimization problem can be expressed as: 

Minimize f(t) = I A.P X. + I A.p.X. (2.l) 

i=1 ^ ^ ^ 

subject to 

- 5 i 0 

tip tip 

u 

o - Of > 0 
root root " 


to 

tip txp “ 


(2.2) 


o ^0 

b.ioot b.root ~ 


u 

b, tip b.tip — 


Z T 1 

m < to < (0^ Wn 

1 — 1 ~ 1 r ^ 


^ C bO x ^ ^ 




Mp - M ^>0 


(2.5) 


and 


Z u. 

X. < X. < X 
3 “ D ~ 0 


If f 


,n 


(2-4) 


where, f(x) represents the weight of the structure and 6 ,a , 
to, Mj^ represent the behaviour quantities maximum deflection, 
maxinum induced stress, buckling stress, natural frequency and 
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flutter Maoh rumber. Eg,. (2.4) represent the geometrical or side 
constraints, which in^sose limits on the size of the design 
variables It can be seen that the objective function, 

Eg. (2.1) as well as the behaviour constraints Bgs. (2.2) are norv. 
linear functions of the design variables. [Ehe constraint, Eq. (2.3 ) 
is not included when the final optimization problan is run. !Ihe 
solution procedure of this optimization problem is discussed in 
Chapter 4 . 



CHAPTER 3 


IIEALIZATIOF OF WIJFG STRUCTURE 


3.1 Stauic Analysis 

The prediction of static, (^namic and aeroelastic behaviour 
of a wing structure represents a problem of extreme complexity. 

The geometric and load conditions are non unifoim, Elementazy 
theories are often incapable of providing accurate results.. Extensive 
efforts have been expended over the years in the development of 
techniques to provide results of acceptable accuracy, and this 
resulted in the development of finite element method as one of the 
most convenient ways to handle complex structures. 

The finite element method of structural analysis, used in the 
present investigation, consists of representing wing structure as 
an assembly of triangular membrane cover elements (9 degrees of 
freedcm), rectangular shear web elements (l2 degrees of freedom) and 
pin- jointed flange elements (6 degrees of freedom) as shown in 
Fig. 3*1» Triangular membrane cover elements carry direct stresses. 

The internal members (spars & ribs) are effective in shear. The 
pin- jointed flange elements represent the direct stress carrying 
capacity of the stringers. The edge displacements are assumed 
linear in all the elements and shear stress constant in the web 
element. The available computer core (32K:, IM 7044 ) eliminated 
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the possibility of using more accurate plate bending elements 
(non- conforming triangular element with 18 degrees of freedom, 
eonfoiming triangular element with 36 degrees of freedom) in the 
place of triangular membrane cover elements. 

23 

Studies made by Gallagher, et al. indicated that the static 
behaviour of a wing structure is predominantly influenced by the 
choice of web element idealization, and not by the choice of cover 
panel idealization. Hence, a comparative study of the displacement 
and eigenvalue results is made to find the rates of convergence 
and accuracies available with different web idealizations. Before 
presenting the numerical results of the comparative study, the 
derivation of the elaaent stiffness and mass matrices will be out- 
lined briefly. 

3.1.1 Derivation of element stiffness and mass matrices 

33ae element stiffness and mass matrices are derived from the 
local assumed displacement or stress field. let the true displace- 
ment state of any element u (x,y,z) be related to the nodal 
displacements U by 

u(x,y,z) = [a(x,y,z)] U (3.l) 

If the strain- displacanent relation is denoted by 

e(x,y,z) = [b(x,y,z)] ^ (3.2) 

iiae matrix [b(x,y, z)] can be obtained by the differentiation 
of the matrix {a(x,y,.z)] . Ihe kinetic energy and potential 
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energy of tbe element are, respectively, given by 

„(e) 1 f ■^3! -v — / \ 

^ ~ J_ P u u dV (3.5) 

V 

and dV (3.4) 

V 

where p is the mass densily, V is the volume of the element 
4 

and n is the velocity. Hie stress vector is related to the 
strain vector e by the Hooke's law 

[ 5 l>[c]t (3.5) 

By assumirig the generalized displacements to be time dependent, 
the potential and kinetic energies of the element can be e^cpressed 
as 

=jU [k] U (3.6) 

and = I U [m] U- (3.7) 

v/here [Jc ] and [ m J are, respectively, the stiffness and equivalent 
mass matrices of the element. By substituting Bg.s. (3* l), (3.2) 
and (3»5) into Sqs. (3.3) and (3.4), we obtain ' 

, = 1 ( / [hf [c] [b] m U (3.8) 

and ^ ( / p{a]®[a] dY)u (3.9) 

By comparing Eh* (3.6) with Eq. (3.8) and Bq. (3.7) with Ifc. (3»9), 
one obtains the stiffness and mass matrices of the finite element 

[k ] = / [b]'^ [c] [bj # 

1 


(3.10) 
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and [m] = J P [a]^ fa] # (3. 11 ) 

V 

Olie integrals in Bis. (3.10) and (3.1l) can be evaluated explicitly 
if the assumed displacement or stress distribution within the 
finite elaaent is known. The detailed derivation of the mass and 
stiffness matrices is presented in Appendix A. 

3.1.2, Itoaerical results for multiweb wing structures 

Several wing structures are analyzed for static deflections 
and natural frequencies by using the following idealizations in 
order to arrive at an accurate and simpler finite element idealization 
for the wing, 

( 1 ) The triangular membrane elements and the rectangular membrane 
elements are used for the idealization of ^ins and webs respectively. 

(2) Both the cover plates and the webs are idealized by the 
triangular membrane elements. 

(3) The cover skins are modelled by the triangular membrane 
elements and the webs are idealized by the triangular shear elements. 

(4 ) Hhe cover plates are idealized by the triangular membrane 
elements and the webs are modelled with rectangiilar shear panels. 

The pirujointed bar elements are used to model the stringeisin all 
the above idealizations. 

The model wings ^own in ilg. 3.2 were tested and analyzed 

23 16 

for static deflection by Gall^er, et al, and also by Rao 




iCtie dimensions of these model wings are given in Eef , 25. iHie 
results for the static analysis and the eigenvalue analysis are 
presented in Tables 3*1 ^md 5*2 respectively. These lesults are 
found to be in complete agresaent with those imported in References 
25 and 16. 

It is observed that the idealization (4) gives results which 
are close to the expaidmental values (Ref, 25) even with a relatively 
coarse mesh size in the case of static analysis, Eventhough no 
esperimental frequencies are available for comparison, the idealization 
vdiich gives lower frequency values is generally considered to be more 
accurate. On this basis also, the idealization (4) can be seen to 
be more accurate. Hence in the subsequent analysis this idealization 
is used, 

5.1,5 Ifeflection and stress a n alysis 

In the optimization problem the wing is assumed to carry a 
specific pay load which in turn is assumed to act as equally 
distributed point loads at the nodal points. In the present worte, 
a clamped boundary condition is specified along the root of the 
wing thereby neglecting any influence of the flexibility of the 
fuselage, Ohe number of degrees of freedom involved in the 
static analysis are reduced to one-half by making the plate 
flexural assumptions. This is made possible by assuming the 
wing to be syrametarf-c about its middle plane amd by choosing the 
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node points on the top bottom surfaces of the wing also ^0 be 
E^ymmetrio* •assuming that (a) -the vertical diepiaifements of 

the upper and lower surface points at a given planfoim location are 
equal, and (b ) the implan® displacsnent of these same respSctive 
points ai-e equal and opposite, the number of the displacement 
variables can be reduced by one-.half. 

Ghe matrix fojiaulation (displacement method) of the general 
structural analysis problon results in the equation 


[K] y = p 


(3.12) 


where [ K ] is the master stiffness matrix of the structure# !Phe 
vectors Y and P represents respectively the dii^lacement and 
the load vectors. Bie ass^bly of the master stiffness matrix 
from the corresponding element matrices and the solution of Bj. (5. 12) 
follow the standard procedure of matrix structural analysis. IHie 
stresses induced in the finite elements can be determined ITma the 
known podal displacements T by using the stress-strain and the 
strain- displacement relations of the linear elastieily. Having 
detemined the displacement and stress field. Dhe maximuin deflection 
of the stiucture and the stresses indiced in the eritioal sections 


of the structure can be restricted by placing: boun^ upon them. 


.u 


- 


< O 


tip. ^ tip "" 




< O 

3?O0t-^ 


(3.13) 


a 


Mp. 'tip ~ ^ 
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3*2 Dynamic Analysis 

ffiae order of the flutter problem may be reduced by introducing 
the first few natural modes of the structure as generalized 

-j g 2^ 

coordinates ^ . Thus the linear eigenvalue problem 

[K ] Y = ^ [M ] Y (3.14) 

is to be solved inorder to set up the flutter equations# Moreover 

bounds are to be placed on the natural frequencies of the structure 

so that anjr possible resonance, that mi^t be caused due to mild 

hamonic forcing, will be avoided# Ohe eigen value analysis refers 

« 

to determining the scalar quantities (eigen values) and the 
corre^onding non trivial solutions Y^^ (eigen vectors) for the 
given master matrices [k]& [m]. Tower method is used to obtain 
the eigen values and the associated eigen vectors. 

For any given degrees of freedom the determination of eigen- 
values and eigen vectors is more expensive than a static solution. 

It is, therefore, desirable to limit the degrees of f3?eedom of the 
already "discrete” system so as to make the eigen- solution economical. 

Hence the static condensation reduction technique suggested by 
25 

Birner, et al, is employed in the present work. The degrees 
of freedom associated with the transverse displacements are retained 
by eliminating those corresponding to the iplane displacaaents. 

Qlius the order of the eigen problem is reduced to one-third of the 
correponding statics solution. 
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If the load vector (Pg) corresponding to the inplane degrees 
of freedom (Tgj is 0 , then the expressions for the reduced stif£. 
ness and mass matrices can be given respectively as, 

= [K^g] 

and [M^] = [M^^]- [M^g] [Kg^l’*'' [K21] - t^22^‘’''^2l^ 

+ [k^ 2^ [Mgg] [K22^"’’ [K^^] (3.16) 

In the case of the reduced stiffness matrix [K^] , Bg,. (3*15), 
none of the structural complexity is lost since all elements of 
the original stiffness matrix contribute. However, in the reduced 
mass matrix, combinations of stiffness and mass elements appear. 

Ihe result is that the original eigen value problem is closely but 

16 

not exactly preserved. Hie natural vibrations of a box beam shown 
in iig, 3*3 are studied in order to see the difference between the 
results given by the reduced eigen problem and original eigen 
problem, Htie results (presented in Table 3 * 3 ) indicate, that the 
first few (s) frequencies and mode shapes obtained from the reduced 
eigen value problaE, are found to be in good agreement with Ihose 
given by the original problem* !Ihis is not considered to be a 
serious problem since in the present work only the first three or 
four frequencies are considered in the flutter analysis. 

3*3 Stability Analysis 

A wing structure may fail in any of the failure modes of 
instability such as over all buckling, skin buckling and stiffener 
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bucfcliug* G-ross buckling of a structure can be analyzed by using 
the concept of geometrical stiffness * . However panel ( stiffened 

plate enclosed between the spars & ribs) buckling is more critical 
compared to gross buckling of wing. Hence gross buckling of the 
wing is net considered in the present analysis. 

16 

In the chapter, a method of including the ^in and stiffener 
buckling failure modes in the optimum design problem is considered. 
Since the covers of a wing structure perform a contouring as well 
as a load- carrying function, the local instability modes of failure 
play an important role in the supersonic regime. Any change in 
the airfoil contour caused by local buckling might lead to an 
inadmissible rise in the drag at high speeds. 

In the modern high performance aircraft the cover skins of the 
wing are usually provided with integral grid type stiffners to achieve 
higher buckling strength for the same solidity. For the purpose of 
present analysis, the portion of the grid- stiffened plate enclosed 
between tihe spars and ribs -is considered to be rectangular and sisply 
supported along the foiir edges, Ihe rectangular plate is assumed 
to have a 0®- 90° grid configuration with longitudinal and transverse 
stiffners of identical shape and spacing as shown in KLg. 3.4, 

3.3.1 Buckling of a stiffened plate 

When the stiffener spacings for the plate are small relative 
to the over all length and width of the plate,- it is possible to 
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use orthotropic plate theoiy to determine the buckling stressi Die 

governing differential equation for the buckling of a simply supported 

27 

flat orthotropic plate under axial compression is given by 


-4 4 4 2 

-n 3 j. -n ^ ^ j. o-n 9 ^ 3 w „ 


(3.17) 


where 


1— V V 


average flexural rigidity of a unit transverse 


cross-section 


T) _ 

2 


aveirage flexural rigidity of a unit 


longitudinal cross section. 




2((xl) = average unit torsional rigidity 


= Poisson's ratio in the directions x and y. 

^ X 


w = transverse displacement 


U = compressive force per unit length. 

A. 


Assuming that plate buckles into one half sine wave and substitution 


w(x,y) = A sin ~ sin 


(3.I8) 


in Eq. (3*17), one obtains 


a 

gross 


' ^ ^ “2 ^2 * 

S' , ■ . 


(3.19) 


33ie smallest value for critical stress is obtained when. 
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( 3 . 20 ) 


and the val’J-e is 


2-n^ 


gross ^ 2 _ ( 


( 3 . 21 ) 


The flexural ^igidifies for the synmietric grid are given by 
_ E ■ 


V2 


I2(l-v ) 


o 4- + r r 

t r, r, (- — 
w s b t \ ^ 


( 3 . 22 ) 


where E = E^ = El^ 


V= V^= 


r, = bn/b» 

b w s. 


V V*e 

and b^, b^| t^, "b^ are the plate and stiffener dimensions as 
shown in Eig, 3.4- • 

The torsional rigidity of the stiffened plate is usually 
considered to be small compared to D^and and the effective 
thickness of the plate is given by 

t = t (l + r r ) 
s s b t 

The theory based on Eq., (3.17) applies to stiffening systons 
which are symmetrical with respect to the middle surface of the 
skin. The as syme try of the stiffening system causes bending 
to proceed, simultaneously with axial loading, and thus the results 
of a buckling analysis for this resembles those for an initially 
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imperfect columii. However, past studies on the buckling of an 

OQ 

integrally stiffened plate with a one sided stiffener system 
indicate that, although seme small amount of bending is evident 
iumediately upon the application of axial load, the ma.TriTTi»m load 
carrying ability of the plate is essentially the same as the 
critical load obtained from ih. (3»2l). 


3.3.2 local buckling of a stiffened plate 


If the skin is thin compared to the stiffens ss, it is 
possible for the individual panel enclosed between the stiffeners 
to buckle as a square plate while the stiffeners undergo a stable 
inplane deflection. In this case, the buckling stress is given by 


2 

4ir E 


local 


12(l - V ) b 


By considering the simultaneous 
29 

failure modes, Gerard derived the 
parameters r^ and r^ as 


(3.23) 

occurence of the above two 
optimum values of the geometric 


^^b^optimum *^*425 

using these optinun values, the buckling stress of the 
isotropic stiffened plate can be expressed as 


o(u) 

bucikle 


= 0,3101 


E 

(l - V^) 



(3.25) 
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3.4 Flutter 

Elutter is defined as a dynamic instability occuj.j[_Qg ^ ^n 
aircraft in flight at a speed called the flutter speed, where the 
elasticily of the structure plays an essential part. 


In most oases, an adequate erraluation of the flutter condition 

is obtained by considering an infinitesimal perturbation about the 

deformed equilibrium position, and analyzing vibration with espo- 

nential tione dependence, since all other motions can be built up 

therefrom by superposition. Hence, theoretical flutter analysis 

usually consists of assuming in advance that all dependent 

variables are proportional to (mreal), and then finding 
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combinations of and « for which this actually occurs 
!Ihis leads to a complex eigen-value problem where there are two 
characteristic numbers which determine flutter Mach number and 
frequency. 


3.4.1 Bormulation of flutter equations 

Lagrange equations in terms of the generalized coordinates 
are used to derive the flutter equations^^’^^ 


dt 



IL. + iL. 

3Ci H± 




i = 1,2,,. 


s 


(3.26) 


where 2 and V are, respectively, the Icinetic and potential 
energies, is the ith generalized force, and e is the number 
of generalized coordinates considered in the flutter problem. 



28 


Itor a discretized structure^ the kinetic and potential 
energies can be expressed as 

T = (3,27) 

Y = ^ C [K] I (3.28) 

and the ^neralized force vector is given by 

P = [q1 I (3.29) 

where [q] is the aerodynamic matrix. 

By substituting Bis. (3.27) to (3.29) into Bi. (3.26), one 
obtains the equations for steady-state oscillations of the 
system, in a state of neutral stability as 


[- [m] + [K] + [Q]]| =0 ( 3.30) 

where a haamionic time dependence, with a circular frequency 
is assmed for the variables 

( 

Derivation of Mass matrix 

Bor a continuous plate-like structure, the kinetic energy is 
given by 

// t P dx dy (3-31) 

where 

i:(x,y) is the thickness distribution, 
p(x,y) is the mass densily per unit area, and 
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w(x,y,t) is the displacement at ary point, 
ibr a discretized structure 

w(x,y,t) = f ^ a ( 3 . 52 ) 

#iere 

W is a column vector of nodal displacements corresponding to 
ttie variable w, and 

a is a column vector of interpotation functions in x and y. 

further, if the deflection pattern can be described by the super- 
position of the first 's' natural modes 

W(t) = [U ] 5 (t) ( 3 . 33 ) 

where 

tu^] is the element modal matrix pertaining to the 
displaconent w(for i^^ element), and 
5 is the vector of modal participation coefficients. 

One obtains 

w(x,y,t) = f ^ [U^]^ a ( 3 . 34 ) 

(Qaus the kinetic energy of the r element can be expressed as 

^ Pt w w dxdy ( 3 . 35 ) 

where dot indicates differentiation with respect to time. 
Substituting Eg., (3.34) into Bg,. (3.35) we obtain 

= -1 [m^] [U^] C (3.36) 

where 

()m 2 = 3- dxdy = discretized mass of 

the i^ element 


(3.37) 
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The kinetic energy of the assembled structure associated with all 
elements can be expressed as 

1 = [u3 ^ [m] [u] I (3.38) 

where [S] is the r x r total mass matrix of the structure, and 
m is the r X s modal matrix, in which the jth column represents 
the jth normal mode of the structure. 

If all the nodal degrees of freedom r are used in Ihe eigen- 
value analysis, the mass matrix [M] of Bi, (3.30 ) is given by 

[m] = [Tif [1] [ul 

sxs sxrrxrrxs 
Derivation of stiffness matrix 

The potential energy of a continuous plat&.like structure is 
given by 

V = ■^ // K^(x,y) w^(x,y,t) dxdy (3.40) 

where is the stiffness influence coefficient at ary point. 

The potential energy of the i^^ element can be expressed 
similar to Bi. (3.56) as 

i 1 T r ii T . i, r i, ->• / v 

c [u ] ^ [K 3 [U 3 e (3. 41) 

where [k^ 3 is the element stiffness matrix given by 

[K^=| // taK^a^dxdy (3.42) 

The potential energy of the assembled stiucture by considering 
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all the displacement -rariables is given by 


7 = i [n ]®tK] [0 ] I 


(3*45) 


Hence the stiffness matrix of Eq. ( 3 .. 30 ) is given by 


[K] = [U ]^ [ K] [U ] 


(3-44) 


Derivation of Aerodynamic matrix 


Second order piston theory is used ia deriving the aerodynamic 
forces acting on the wing. !Ehe .basic assumption of this theory is 
that the pressure at any point on the wing depends only on the 
normal component of fluid velocity at that point, and is related 
to the later in the same manner as Ihe pressure behind a piston 
moving in a one-dimensional cylinder. Obviously, 3-dimensional 
effects are not included in this theory. However these effects on thin 
wings are relatively snail at hi^ flight speeds. 


According to this theory, the pressure differential is given 


by (with the factor e^**^ suppressed) 


V 

CO 


Y +1 ■« 9Z 


yx,y) = (1 + ^ ^ 


. ) (v w) 

9x ® 3x 


a 

00 


(3.45) 


where 

= mass density of air 
a = speed of sound 

Y = ratio of specific heats = 1,4 
ir = free staream velocity 

22(x,y) =fthic£hess distribution of wing as shown in ilg. 3*5 
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tte aeiodynajic matrix ^ 

thrcagh a aalcalatioa of ririusl woifc. ■ Iha virtaal K,ac for the 

•■til . 

1 element as given by 

= // w Ap dxdy 

= 2 2. 2-.^„ " If ^>=^7 + 2 „(y+ 1 // ; II If arfy 

+ i a P.a» /; ; w drfy + 1 p_(y+i) // w ^ w dray 

( 3 . 46 ) 

where w is the virtual displacement and the integration is over 
the area of the element. 

Introducing Eq. (3.32) into (3.46) yields 

[Q^] i 

where the required aerodynamic matrix (whose order is equal' to the 
number of nodes of the ith element) [Q^] is given by 

[Q^-] = 2 P^a^ v„ [a""] -H2ip^a^« [B^] + 

Poo (y+i) vf [o'"] i p„(y4.i) v„ m ] (3.47) 

pidth 

[a"i = //t" || dfdy 
[b"-] = ff t^S drdy 

[Bh = // ||J"t drfy 


( 3 . 48 ) 
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In the present woik, down stream slope of an element is taken as 
the average slope of its nodes. Also triangular membrane elements 
represent the planform elements and the vector of interpotation 
functions for these elements is given by 


a (x,y) = 


2A, 


123 


(y^- yg) (x - x^) - (x^- Xg) (y - yg) 

^(y^~ y^) (x- x^) - (x^- x^) (y- y^)y (3.49) 

(yg- y^) (x- x^) - (xg- x^) (y- y^^ 


where A^^j area of the eluent shown in Pig. (3.6). Bie 

double integrals of Bq.s. (3.48) are evaluated using the following 
integration formulae. 

If the origin is at the centroid of the triangle, then over 
the area of the triangle 


/ / dxdy = A^23 

/ / xdxdy = 0 

/ / ydxdy = 0 

/ / x^dxdy = + X3) (3.50) 

/ / y^dxdy = (y^ + yg + y^) 

A , 

/ / zydxdy = — (x^y^+ X2y2+ x^y^ ) 


Since only a translation of X, T axes is involved, the trans- 
formation matrix is sin^ily an identity matrix. 
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!Ehe aerodynamic matrix for the assembled structure is obtained 
from an assemblage of the individual plan foim element matrices 
[ 3 by the standard procedures of structural analysis. Hius 

the airforce matrix to be used in Bj. (3.30) can be obtained as 


= 2 2 i a^m [B] + (Y+i) v^ [C ] 


+ i P^C>+1 ) V £0 [l] 


where [ A ] 


Eu] [A^][u3 


(2.51) 

(3.52) 


[D] = [U3 " [D^3 [U] 

5.4.2. Solution of the flutter problan 

IThe requirement for non trivial solution to Bq. (3.30) is that 
the deteiminent of the coefficient matrix of ^ must vanish. Baus 
the flutter equation becomes 


I [Z] - [M] + IQ ] 1 = 0 


(3.53) 


Bie Sq. ( 3 * 53 ) represents a complex, nonlinear, double eigerv- 
value problem since there are two unknowns, free stream velocity 
V and oscillation frequency uj . Bq. (5.53) can be written in a 


more convenient form by noting that 


[M] 


and 


Ez ] = fNyJ 


(3.54) 


when the noimal modes of vibration are used as generalized 
coordinates. Substituting Bqs. (3.54) into Bq. (3.53) and with 



35 


some reairargement, one obtains 


I- (. [A]+^|-tc] ) 


00 00 


1 


with 


+ i (^) ( [B] J- CD] ) I = 0 

a) ^ 


(3-55) 


~ reduced frequency, and 
= some reference length. 

!Ihe Bq. (3-55 ) with Mach number (~) and reduced frequency (k ) 

“w ^ 

30 

as the unknowns is solved by V-g method , a double iterative 
process. 


3.4.5 numerical example . 

In order to test the accuracy of the flutter equations and the 

flutter analysis program, a double wedge airfoil shown in Pig, (3.7) 

is analyzed for flutter. !Ehe supersonic flutter analysis of this 

13 

airfoil was previously reported by Schmit and Thornton . In this 
work authors assumed a beam- type structural behaviour ln determining 
the bending- torsion flutter characteristics of the wing. 

l6 

In the work reported by Rao , the same problem is analyzed 
by some finite element formulation, with 40 triangular membrane 
elements as shown in Pig. (3.6), 55 degrees of freedom in eigen» 
value analysis and 4 degrees of freedom in flutter analysis. In 
the present work the same problem is analysed with 15 degrees of 
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freedom in eigenvaJne analysis and 3 degrees of freedom in 
flutter analysis. !Ihe results are reported in !Eable 3*4. The 
present analysis underestimates the flutter Mach number by about 
5^. This can be considered as good agreement since the degrees 
of freedom in present flutter analysis are only 3, 
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TABLE 3.1 

DISPLACEMENT RESULTS POE ilULTIWEB WING STRUCTURES 
Uniform, total load applied = 1500 lb. 


!]^pe of wing 

Maximum tip 

deflection in inches 


structure 

idealization( 1 ) 

ideali- 

zation(2) 

ideali- 
zation(3 ) 

idealii 
zation(4 ) 

Experiaental 
(Ref. 23) 

Model 3 

0.445 

0.240 

0,280 

0.704 

0.797 

(sweepbaok 45°) 






Model 5 

0.613 

0i303 

0.316 

1.345 

2. 120 

(sweepbaok 30°) 







# 

TABLE 3*2 




EIGENVALUE RESULTS POR MULTIWEB WING STRUCTURES 


Type of wi23g 
structure 

First natural frequency 
(rad, /sec. ) 

Second natural freruenoy 
(rad, /sec. ) 


ideali- ideali- 
zation zation 

. (1.) . .(3). 

ideali- 

zation 

(4) 

ideali- 

zation 

(l) 

ideali- 

zation 

i3)„.„. . 

ideali- 

zation 

(4) 

Model 3 

(sweepbaok 45°) 

578 725 

468 

1482 

1594 

1392 

Model 5 

(sweepbaok 30 °) 

606 ■ 846 

416 

1522 

1676 

1402 




TiBEE 3 


MTDEisjj lEEQUEHClES OP 


Humber of the 
Eigenvalue 


Hatural frequej 


60 degrees 
of freedom 


1 

2 

3 

4 

5 

6 
7 


43-0 

88.5 

170.8 
235.5 

318.9 
395.4 
455.7 


8 


504.3 


9 


602.5 


10 


704.1 
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(rad« 


20 degrees 
of freedom 


mode shape 


43-8 

88.7 

171.2 

235.8 

319.1 

395.7 
456.0 

504.7 


first bending mode 
!Eorsional mode 

second bending mode 


Inplane mode 


704.2 
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Talile 3 4 

HESUISS OE ilEUTTER MilffSIS jPOH A IXDUBLE WEDGE AIREOID 

Attitude = 30,000 ft, material = steel 

T = 1.50 ft, C = 7.50 ft, t = .058 ft 

c 


Naturel frequencies Flutter frequency 
(rps) (cps) 


Flutter Mach number 

Ref .13 Ref .16 Present analysis 


44.60 

201 .46 

217.53 


22.65 11.77 13.57 14.10 

-11.99 +545 

io error io error 


219.06 





MOCEL 3 MODEL 5 


FiG-3-2 MULTIWEB BOX TYPE WING STRUCTURES 

( Geometry and node point locations and 
other dimensions in reference 23) 





Mate rial -Steel ’ ■' 

Cover thickness: 0-0.5 

Web thickness-. Q-Of ■ ■' 

Idealization: cover plates - triangular^;g 
membrane elements'^- 

Webs - shear triangles 


^ 




FIG-3;3' A:;^ 


BEAM 








Non- dimensional' parometers, 'n 


FIG-3'A A STIFFENED PLATE UNDER AXIAL 
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CHAPTER 4 


OPTDIIZATIOir AJLGORITHM 

The problan of optimum design of aircraft wing structure has 
been foiiaulated. as a nonlinear mathematical programming problem# 

The selection of solution procedure and its description is presented 
in this chapter. 

The solution procedures of a math an atical programming problem 
can be classified into (a) Direct methods, in v/hich the constraints 
axe explicifly handled, and (b) Indirect methods, in which the 
constrained foimiilation is transformed into a sequence of unconstr- 
ained optimizations. One main reason for the appeal of the 
sequential unconstrained optimization foimulation of the constrained 
optimization problem is that the sequential nature of the method 
allows a gradual approach to the criticality of the constraints. 

In addition the sequential approach permits a graded appiiszimat ion 
to be used in the analysis of the system. Also this transformation 
avoids the necessity of coping separately with the boundary of the 
inequality constrained region, e#g., by attempting to move along 
the boundary once it is encountered, Shch a motion is cxmibersome 
when the constrained surface is nonlinear and it requires the 
solution of anothar prograaming problem to find feasible direction. 
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One of the approaches to reduce the total computational time 
of automated optimization problems is to adopt a method which 
peimits the use of approximate analysis without involving signi- 
ficant errors. Ehe penalty function methods allow the use of 
approximate analysis during the early phases of the optimization. 
Purthermore the variable metric unconstrained minimization 
technique, discussed in a later section, is inherently more stable 
and little effected by minor errors introduced through analysis 
approximations. It is for this reason that an optimization 
technique from the indirect methods was choosen for the purpose of 
the present work, 

4.1 3?iaoco.-McGoimick Interior Penalty lUnction Method 

Penally function methods transfoom the basic probloE into 
alternative formulations such that numerical solutions are sought 
by solving a sequence of unconstrained minimization problems; A 
penalty function ^(X , r^^.), will be formulated by adding a 
penalty term *)» which is a function of the constraints 

to the objective function f(x) and this penalty function will be 
minimized for a sequence of response factors r^^. There are many 

^ . 35’ 

choices available for the penalty term, however Piacco-McCoimick 

suggested a simple form rjj. I In this formulation the 

j=1 

penalty term is small at points away from the constraints in the 
feasible region, but it "blows up" as the constraints are approached. 



Taen the transformed optimization problem is 
Minimize 

$(x.r^) = f(J) + I (4.1) 

3=1 gj-(x) 

for a decreasing sequence of r, , r < r . 

k' k+i K 

Tae penalty tem is not sin^jly defined if x is infeasible. In 
mar?r engineering problems it is not difficult to find an initial 
feasible point at the expense of large value of f(x). If there 
is ax3y difficulty in finding an initial value of X, the method 
described in reference 34 could be used to fird a feasible starting 
point. In the present woifc the feasible starting points are 
found by a process of trial and error. 

Since each of the designs generated by this method lies 
inside the aceeptible region of the design space, the method is 
classified as an "interior penalty function" formulation. She 
method tends to generate a sequence of designs vshich decrease the 
value of the objective function such that none of the designs in 
the sequence is critical with respect to the set of inequality 
constraints. aJhis characteristic makes it possible to use 
approximate analysis methods during major portions of the optimi- 

iq 

zation procedure 

4.2 Bavidon-gletcher^Powell Variable Metric Unconstialned 
Miniinization Method. 

The selection of an efficient unconstiained minimization 
algorittni is extranely iaportant since a sequence of such 
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miniiaizations of (?> (x, r^) are to be performed. Biere is wide 
choice of algorithms of search. These methods are broadly 
classified into (a) methods of direct search vshich use only 
function values and (b) gradient methods which use the gradients 
of the fmction also. In general the gradient methods ai« more 
superior to the non-gradient methods since they use more infor- 
mation about the function (namely, the gradients). Tbr this reason 

36 

the DavidoruZLetcher-Powell variable metric method is used in the 
present study. This algorithm is probably the most powerful 
general procedure mvi known for finding a local unconstrained 
minimum of a function of mary variables. 


* 4 * 

Given a starting point and a possitive definite matrix 

[h 3 , this method seeks the local minimum of <^(x, r^^) by 
0 


generating a sequence of vectors i=1,2,,., 


where 


^i 


and = - [H^J V (x^) 


( 4 * 2 ) 

(4.3) 


Thus X is the design vector corresponding to the minimum of 
i+1 


<{>— function along the current direction S^, > ihe starting 

design vector and t* is the minimizing- step length. The i^ ^ 
iteration begins with the deteimination of the search vector 
according to the relation (4«3)# 

Prior to begining another cycle in the iteration, the matrix 
iH-l is modified to taie into account local characteristics of 



50 


"the $ -•fliGc'fcion inorder "to avoid the zig-zag behaviour cominoii to 

other mi n i ini zatioti teohniciues. !Ehe procedure for the 

updating of the matrix [H^] is given in reference 36 . Ihe 

stability of the procedure is insured by preserving the symmetry 

and the possitive definiteness of [H^] viiile it is updated. Bie 

positive definiteness of the [H^] matrix is influenced only by 

the accuracy with which t* is determined. In applying the 

algorithm, therefore, care must be taken to insure that the [H^^] 

matrix is not updated with data arising from poor approximation 

^ 0 } , 

to T*. Eierefore, whenever la3?ge ( x* is poorly 

approximated) the on^dimensional minimization algorithm may be 
reapplied to refine x*. If tiiis procedure takes excessive 
computational effort, (more than 2 or 3 refits) the updating process 
is skipped ( set [ ] ) and the algorithm is continued 

as before. 

4.3 Cubic I nt erpotation Ona-Pimensional Search Technique . 

Several methods are available for determining the minimization 

step lengths t* in (4.2). Since most of the effort goes 

towards the computation of minimization step lengths, the Belection 

of an efficient one-dimensional minimization technique becomes 

extremely ia^iortant* In efficient one-dimensional search method 

which is a natural, choice with the variable metric algorLiisn. is 

34 

the well known cubic inteipolatLon one-dimensional technique . 




51 


The cubxc xateirpolation technique is essentially a gradient 

algoritha and the gradients already computed in the variable metric 
algoritha can readily be used. 

She one-dimensional search problem is to find t = i* which 
yields the first local mininum of 


+ S^t) = <j,(T) 


subject to 0 < T < 


(4.4) 


where t is the largest t for which t ) lies within 

the feasible region. 

She minimizing step length is given by 


$4 + Q _ z 

* 7 . -D 

T =Us 






(4.5) 


where 


4(t=. t^) = 4 ,^ 

4 .(t= Tjj) = 


d4> 

ST 


T = TA 


(4-6) 


M- 

9t 


T= T, 


= (^’ 
B 


and (j) • < 0 < $ ’ 

■ A ' ''.S 

also z = - — + 4,^ + ^4 


T — T 

B A 


Q =. [Z^- 4^ 4;^ ] 
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Ihe effort involyed in this algorithm is reaching the points 
T = A, aM T = B satisfying the condition 
initial step length is small numerous increases in x -will be 
necessary before this condition is satisfied. On the other hand 
if the initial step length is chooser comparatively large, then 
the interpolation takes place over so large an interval that it 
produces a poor approximation to x*, A sinple solution to this 
problem is to use a priori method which assumes initially that 
(})(x) can be approximated by a quadratic and then use 4>(o), <!>’(o)» 

, a guess at the minimum a£ ^ along S as the data for 

interpolation. For <{> to be minimum, the initial step length 

34 

can easily be calculated as - „ . 

5)A’ U-'?) 

In the present work 5 is assumed to be £5^ of $(o), this 

being an educated guess# 

nhe one-dimensional interpolation procedure is terminated 
when the cosine of the angle between and is anall, 

i±l <e(=0.05) ( 4 - 8 ) 

If this orthogonaliV test fails, the interpolation is again 
performed setting 
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(4.9) 


(4.10) 


In the present -wooi: these refits are limited to a maximum of 5 in 
each one dimensional search problem, 

4.4 Additional Gonsideratioos 

( 1 ) Starting Point 

!EhG feasible starting point is found by a process of trial and 

error. Each subsequent stage uses the solution of previous stage 

as a starting point. In some cases, the overall procedure has been 

34 

accelerated by employing an extrapolation technique to determine 
startling points for subsequent unconstrained minimizations. Starting 
points obtained by extrapolation must be checked to be sure that 
they satisfy the constraints. 

( 2 ) Values of 

If r is very large the function is easy to minimize, but 
the minimum may lie far from the desired solution to the original 
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constrained minimization problem. On the otherhand if r is 
small the function will be hard to minimize. In the present work, 
the value of r^ is choosen such that 

1.25f (X^) < 2.00 ( 4 , 11 ) 

3}he subsequent values are found by using the ratio 

( 3 ) leimination of minimization for each r . 

K 

The minimization of <!*(x,rj^) for each 3 ^ is terminated when 
the deca?ease in penalty function for successive iterations is less 
than 0.5/2. To gua 3 ?d against premature termination due to slow 
convergence a mini m um of IT (no. of variables) one dimensional 
minimizations is performed before testing for convergence. Even 
after N iterations if the convergence is failed, the [h] matrix 
is reset to [ 1 ] matrix and the minimization process is continued. 

(4 ) Relative minima. 

In order to see whether aiy relative minima exist in the 
design space, two completely different starting points are used^^^^^^ 
for the sequence of minimizations for one problem. The 2 sequences 
led to the same optimum design (except for a ^all difference that 
might have occured due to numerical instability) nnd hence it 
appears that the local optinum is the same as the global optimum 
for that problem , 
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(5) Reducing the total computational time. 

It has been obserred that some of the automated optimum design 
problems ucfee a longer time to satisfy the prescribed convergence 
criteria even after reaching veiy near to the optimum design. !Ehis 
happens whenever the function is highly distorted or eccentric. 

In such cases it is HKJt worthwhile to try to reach the exact minimum 
to obtain about 1 or 2^ decrease in the obj ective at the expense 
of 40 to 50fo more computing time. This problem can be tackled by 
having coarse convergence criteria in the early stages of optamiza- 
tion (initial 'r' s) and later using refined convergence criteria. 

Also initial unconstrained minimizations produce greater reduction 
in the objective function and the final unconstrained minimizations 
produce lesser reduction in the objective function as the design 
point gets closer to the boundary of the constraints . Thus it 
may not be profitable to press for a few percent reduction in the 
later stages of unconstrained minimizations at the cost of substantial 
computer time (30 to 40^). Hence in the present work only 4 
unconstrained minimization are employed for each problem. In 
such a case testing for the satisfaction of Kuhn- Tucker conditions 
can be deleated. 

Also most of the computer time in the automated optimum design 
programs is expended in repeat tini®~consumming design analyses. 

Hence a rapid reanalysis is the key factor in speeding the optimi- 
zation program. The behaviour variables can be linearly approximated 



whenever the changes in the design variables are small. 

let Mj, be the flutter Mach number corresponding to the design 

* 

X and. Mp be the flutter Mach number corresponding to a perturbed 
design X . % knowing M^, and the rates of changes of Mj, at 

• 4 * 

X y the flutter Mach number at the perturbed design can be 
approxi m ated as 

N 3M 


* 






( 4 . 13 ) 


urovided 




< e , k e U 


( 4 . 14 ) 


!rhe value of e is taken as 0.1 in the early stages of optimization 
and is successively reduced to 0.02 in the final stages of optimization 
A similar procedure is ^plied to compute the other response variables 
like natural frequencies, deflections & stresses. It has been 
observed that this si^le approxi m ation alone has resulted in 30^ 
reduction of the total computer time, however with about 2^ increase 
in the final minimum we ight , 

(6) Partial derivatives, 

3he partial derivatives of the response variables, weight, 
natuaral frequencies, deflection, flutter Mach ramiber are computed 
using foirward finite difference formula, for example. 
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where M«. is the flutter Mach number corresponding to the design 
vector X* = 1 , 05 ^ , 

Although exact expressions for the partial derivatives of the 
behaviour quantities are available , their usage requires additional 
storage in the computer programme and hence the partial derivatives 
are computed using the finite differs me method. 



CHAPTER 5 


RESULTS Aim DISCUSSIORS 

A computer programme packagehas been developed to obtain the 
minimum weight design of a symmetric wing (airfoil section is of 
diamond shape) with constant planform area satisfying the require- 
ments of strength, stability, frequency and flutter (not included in 
subsequent study), !I!he programme is written in Eortran 17 language 
and epcecuted on IBM 7044 computer system. Also a flow chart of the 
same programme is presented in Appendix B, lihe dimensions of the 
wing structure (Pig. 5.l) depend upon the aspect ratio, thickness to 
chord ratio and sweepback angle which are among the design variables 
considered. The leading edge and trailing edge spars are assumed 
to be at 10^ and 90^ of the chord respectively. Then the main 
dimensions of the wing can be obtained as 

b, semxspan = (AR . SA , 0,5)”*^^ (5.l) 

c, chord = b Tan (SW) + SA/b 

t, max imam thickness = (c - 2 b Tan (SW)). ffi./0,8 

vdiere 

SA — Ilanfom ai*ea of the semispan of the wing, 

Sy= Sweepback an^e of the leading edge, 

AR « Aspect ratio of the wing, 

OK =; Thickness to chord ratio of the wing, 

Bie relevent data assumed for the wing is presented in table 5.1, 
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As sbown above, the main dimensions of the wing depend upon design 
variables and whenever the design is changed, the dimensions are to 
be re-evaluated. 


ISae finite element idealization of the wing is done with 50 nodes 
and 172 elaaents (64 triangular membrane elements, 36 rectangular 
shear panels and 72 pin-jointed bars). !Che original 120 degrees of 
freedcmt are reduced to 60 using plate flesxural assimptions, and these 
60 degrees of freedom are further reduced to 20 in the eigen problon 
by making use of static condensation technique. The flutter 
constraint is not included in the optimization problem since its 
double eigenvalue analysis is time consuming and also it requires 
additional storage in the computer system. To start vri-th 13 design 
variables - 10 design variables representing the linearly decreasing 
mass distribution of cover skins, spar webs, rib webs, spar flanges 
and rib flanges, and 3 more design variables representing the aspect 
ratio, thickness to chord ratio and sweepback angle were considered. 
However the results of the initial paremetric study (Table 5.2) 
conducted to observe the influence of the design variables on the 
behaviour quantities showed that the linearly decreasing mass distri- 


bution of spa 3 ?s and ribs can be r^laced by uniform distribution 
without any appreciable change in behaviour constraints, thus reducing 
the number of design variables to 9. and X 2 represent the 
linearly decjr easing mass distribution (X 2 *" xX^) of the cover plates; 
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and Xg represent the flange areas of spars and ribs and 

^ represent respectively the aspect ratio, thickness to 
chord ratio and sweepbadc angle. 


!Ifhe minimum weight design problem is initially solved employ- 
ing linearly - approximated redesigns and the results are given 
in Sable 5.3a. So assess the saving achieved in computational 
effort by employing the linearly - ^proximated redesigns, the 
same problem is solved (with the same starting point) using exact 
analysis in the third and fourth unconstrained minimizations. The 
restilts (o&ble 5.3b) indicate that the total computational time is 

reduced by about 30^ and the minimum weight is increased by about 

. %CO-\<X^hhCf> 

1,50% when the linearly - approximated redesigns are used. 


To check whether the optimum obtained is local or global the 
same problm is solved with a different starting point. The 
results are presented in Gfeble 5.3c, It is observed that this 
second starting point gave an 8% over-wei^t design, suggesting 
another local optimum. In such a case it is customaiy to solve 
the problem from another starting point and obtain some inference 
regarding the global optimum from the three local optima by using 
inteiTpolation techni<iue. However this is not attempted in the 
present investigation. 

It can be observed that the optimum structure, in the present 
case has become stiffen, Ihis can be justified with the argument 
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that the present idealization reduces it to a Mitchell structure 
aiai an optiirm Mitchell structui*e is a stiffen stiuoture. 

aSae results of the parametric study conducted about the 
optimuai point of this probloa are presented in Mgs, 5,2 to 5.9. 
!Hiese graphs help to evaluate the extent by which one pays the 
penaligr interns of wei^t if off - design points are considered. 

It can be observed from the results ( (Table 5.3) that only the 
root buckling stress constraint, and no other side constraint, is 
active, Normally it is expected that an optimum solution should 
be a lower bound solution. However, in the present - case a wide 
margin is given on the side constraints, which explains why none 
of them are active, (Do check whether the optimization routine is 
working correctly or not, a second optimization problem is solved 
with constraints made more stringent. Hie results of this problem 
are presented in (teble 5.4. The optimum of this problem is 
characterized by an active frequency constraint and some side 
constraints. 

]i\irther, it can be observed from the results, that minimum 
wei^it design of a wing structure is feasible even when the 
geometric parameters aspect ratio, thickness to chord ratio and 
sweepbaok angle are included in the design variables. 
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DATA 

TABLE 5.1 

FOR EXiiI>!PIE . F/ING 


Material properties: 

Material 

Titanium 


Young's modulus 

6 

16.4 X 10 psi 


Poisson's ratio 

0.3 


Density 

0. 16 1b./cu,in 

Details of wing loading: 

Planfoim area 

3542 sg. ft. 


Pay load 

400,000 lb 


Solidity ratio 

0,1 (ist problem) 



0.11 (2nd problem) 

Flight conditions: 

Alititude 

25,000 ft. 


Density of air 

2 ^ 

0.0010663 lb - sec /ft 


Sneed of sound 

1016 . to ft/sec. 


Flight Mach number 


2.5 


RESULTS OP IIITIAL PARAMETRIC STUDY 






TABLE 5.2 (contd 
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I ABLE 5.3a 


OPTMIZATIOE EESUITS POA TrIE FIRST PROBIEM 
(Employs linearly approximated redesigns whenever AS^<0.1) 


Behaviour 

Quantities 

Lower 

Bounds 

Upper Initial 
Bounds Design 

Sequential Unconstrained Minimir,- 
atiors 





r ^=2 . 9 - 

r 2=0 . 2 9 

r ^= 0 . 02 9 r 

6 ( in) 

— 

69.0 

46 .0 

37.6 

36.34 

36.34 


- 

75.0 

42 .0 

43.8 

49.31 

49.31 

‘’tip 

- 

75.0 

8.03 

8 .66 

9.82 

9.82 


~ 

55.0 

41.9 

43.7 

49.19 

49.19 


- 

55.0 

7.31 

7.97 

8.30 

8.30 

0)^ (cps) 

1.2 

4.5 

2.30 

2.96 

3.45 

3.45 

“2 (®p® ) 

2.0 

7.5 

3.80 

4.91 

5.23 

5.23 

Penalty Function 


103500 

87000 

34500 

29500 

Weight (lb 

.) 


51900 

38300 

29100 

28990 

lime 



about 

40 : minutes 



* denotes active constraints 
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TABLE 5.3a (contd.) 


Design 

lower 

Upper 

Initial Sequential Unconstrained Minimizations 

Variab- 

les 

' Bounds Bounds 

Design 

r^=2.9 

r2=0.29 

r_=0.029 

3 


0.01 

0.06 

0.035 

0.0259 

0.0193 

0.0193 

^2 

0.10 

0.60 

0.35 

0.2590 

0.1925 

0.1915 

X, 

0.10 

0.40 

0.20 

0.1988 

0.1750 

0.1750 

^4 

0.10 

0.50 

0.35 

0.1960 

0.2790 

0.2785 

Xc- 

5 

0.10 

0.40 

0.20 

0.1988 

0.1750 

0.1750 

^6 

0.10 

O 

• 

o 

0.35 

0.1960 

0.2790 

0.2790 

X? 

0.75 

3.0 

1.50 

1.491 

1.3120 

1.3120 

^8 

0.01 

0.05 

0.015 

0.0189 

0.0190 

0.0190 


10.0 

35.0 

20.0 

21.40 

21.50 

21.50 

Number of one 

minimization 

unconstrained 

-dimensional 
for each 
minimization 

9 

4 

2 
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TABLE 5.3 b 

OPTMIZATIOF RESULTS EOS THE FIRST PROBLEM 

(Employs exact re-analyses in 3nd and. 4tli uncostrained 

minimizations) 


Behaviour 

Quantities 

Lower 

Bounds 

Upper 

Bounds 

Initial ' 
Be sign 

Sequ-ential Unconstrained 
Minimizations 

r^=2.9 r2=0.29 r^=0.029 

r'^=0.00 
4 29 

6 (in) 

— 

69.0 

46.0 

37.6 

36.34 

41.7 

41.7 


- 

75.0 

42.0 

43.8 

49.31 

54.4 

54.4 • 


- 

75.0 

8.03 

8.66 

9.82 

10.5 

10.5 


- 

55.0 

41.9 

43.7 

49.19 

54 . 2 

54 .2^- 



55.0 

7.31 

7.97 

8.30 

8.96 

8.96 

(cp^) 

^•2 

4.5 

2.30 

2.96 

3.45 

3.24 

3.24 

“2 - (CP3.) 

2-^0 

7.5 

3.80 

4.91 

5.23 

5.03 

5.03 


Penalty 

Funct ion 

103500 

87000 

34500 

29100 

28380 

Weight (lb.) 

51900 

38300 

29100 

28300 

28300 


Total Time about 60 , minutes 
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TABIE 5.3b (contd.) 


Des j^n 

Lower 

Ujoper 

Initial Sequential Unconstrained Minimations 

Variab- 

les 

• Bounds 

Bounds 

Design 

ri = 2.9 

- r 2 = 0 . 29 

jl 

0 

• 

0 

9 r ^= 0.0029 

Xi 

0.01 

0.06 

0.035 

0 , 025 S 

0.0193 

■ 0.0188 

0 . 0188 


0.10 

0,60 

0.35 

0.2590 

0,1925 

0.1880 

0.1880 

X3 

0.10 

0.40 

0.2 

0.1988 

0,1750 

0.1818 

0,1818 

-4 

0.10 

0.50 

0.35 

0.1960 

0.2990 

0.2780 

0.2780 

-s 

0.10 

0.40 

0.20 

0.1988 

0.1750 

0.1820 

0.1820 

■^6 

0.10 

0.50 

0.35 

0.1960 

0.2790 

0.2780 

0.2780 


0.75 

3.0 

1.50 

1.491 

1.3120 

1.368 

1 . 368 

Xq 

0.01 

0.05 

0.015 

0.0189 

0.0190 

0.0191 

0.0191 

Xe 

y 

10.0 

35.0 

20.0 

21.40 

21.50 

20.20 

20.20 

Ifiiiaber of one- 

minimizat ions 

unconstra,ined 

•dimensional 

for each 

optimization 

9 

4 

8 

2 
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TABLE 5.3c 

OPTIMIZATIOH RESULTS FOR /THE FIRST PROBLEM 
(starting point is different) 


Behaviour Low- 
Quantities er 

Boun- 

ds 

Upp- 

er 

- Boun 
ds 

Init- 

ial 

-Des- , 
•ign ' 

Sequential Unconstrained Minimizations 

^1= 

BOO 

: r2= 

80 

r 

8 

3“ 

^4“ 
0 6 

^5T 

O.Od < 

^6~ 

0.008 

6 

( in) _ 

69 

.0 

39 

.00 

41. 

,08 

38, 

.80 

39. 

.64 

40. 

29 

42, 

.85 

42. 

95 

‘^Root 

(ksi) - 

75 

.0 

45 

.00 

47. 

.01 

45. 

.93 

48. 

.29 

50. 

71 

52. 

.41 

52. 

.47 

^t ip 

(ksi) — 

75 

.0 

8 

.88 

9. 

.07 

8, 

.75 

8. 

.87 

8. 

79 

10 

.13 

10, 

.19 

'^b .Root 

,(ksi) - 

55 

.0 

44 

.30 

46, 

.87 

45. 

.80 

48. 

.14 

50. 

55 

52 

.25 

52, 

.32 

*^b . t ip 

(ksi) - 

55 

.0 

7 

.45 

7. 

.72 

7 

.34C 

1 7 

.62 

7. 

56 

8 

.75 

8, 

.81 


( cps)l.2 

4 

.5 

2 

•92. 

2 

.89 

2 

.95 

2 

.97 

2. 

.95 

3 

.07 

3 

.08 

“2 

( cps) 2.0 

7 

.5 

4 

.48 

4 

.49 

4 

.58 

4 

.71 

4. 

.74 

4 

.83 

4 

.84 


/ 

Penalty Function (xlO^) 12.5 1.275 .158 .046 .0324 .0309 

Weight (lb.) (xlO^) 38.0 36.5 36.0 34.7 33.4 31.2 30v5 


Total Time 


about 90 minutes 
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TABLE 5o^o (contd.) 


Design 

Lower 

Upp- 

Init i“ 

■ Sequential U 

nconstrained 

Minimi 

.zations 

VariaL- 

les 

- Bounds 

er al 

Bounds Desi- 
gn 

CO 4 

0 H 
0 11 

r2= r 

80 8 

3 = ^ 4 = ^ 

.8 . 

5= "*6 = 

08 .008 


0.01 

0 e 06 

0,042 

,029. 

... . 028 . 

.024 

.i019 

.025 

.025 

X2 

0,10 

0.60 

0.320 

.253 

.251 

.235 

.220 

.215 

.213 

x. 

0.10 

0.40 

0.220 

.202 

.200 

.192 

.182 

.183 

.184 

""a 

0.10 

0.50 

0.210 

.213 

.210 

.215 

.214 

.212 

.213 

=^5 

0.10 

0,40 

0.180 

.187 

.186 

.190 

.192 

.190 

.190 


0.10 

0.50 

0.368 

.209 

.191 

.214 

. .205 

.198 

.199 

B’ 

0.75 

5.00 

1.550 

1.410 

1.595 

1.430 1.440 

1.430 

.l-,43 

^■8 

0,01 

0.05 

0.017 

.018 

.0183 

.:0189 .0195 

.0193 

.0192 

^9 

10.00 

35.00 

16.000 

17.8 

17.5 

17.8 

17.9 

17.8 

17.8 


number of one- 

-d ime ns io nal 9 

7 

11 

4 

4 

2 

minimizat ions 

for each. 






unconstrained 

optimization.. 
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TABLE 5-4 a 


OPTMIZATIOI ILSSULTS FOR THE SECOID ■ PROBLEM 


Behavioiir Lower Upper Initial Sequential Unconstrained Minimizations 
Quant it i- Bounds Boun-Design r-=0.95 r^=0.095r ^=0.0095 r. =0.000 
es (is J- ^ p 4 g 5 


6 ( in) 

■’ 47.0 

45.96 

40.14 


41.37 

41.57 


“ 75.0 

42.05 

45.66 

CO 

46.78 

46.78 

®Iip 

- 75.0 

8.029 

8.968 

03 

CD 

?H Jh 

o 

9.111 

9.111 


■“ 66.5 

41.91 

45.54 

fXU 

CO 

si -H 

46 .66 

46 .66 


).Tio^^®^^ •* 66.5 

(cps)l.OO 3.00 

(cps)3.00 5.00 

7.311 

2.302 

3.794 

7.214 

2.965 

4.282 

Optimizat ii 
down for t; 

7.334 

2.93 

4.250 

7.334 

2.93^^ 

4.250 

Penalty Function 

102000 

79000 


37000 

37000 

Weight (lb. ) 

51000 

37400 


36960 

36900 


Total Time 


About 60 minutes 
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TABLE 5.4a (contd.) 


Lesign 

Lower 

Upper 

Initial Sequential Unconstrained 

Minimizat io ns 

Variab- 

■ Bounds Bounds 

Design 

r^=0 .95 

r2=0 

.095 r^=.00 

95 r^=0.0g|- 


0.02 

0.04 

0,035 

0.0253 

- 3 

0.0250 

0.0250 

^2 

0.20 

0.40 

0.35 

0.2530 

XT 

0.250 

0.250 

^5 

0.15 

0.25 

0.20 

0.1970 

1 

0.1932 

0.1932 


0.25 

0.40 

0.35 

0.330 

CQ 

05 

0.328 

0.328* 


0.15 

0.25 

0.20 

0.1970 

0 

Fh 

rO 

0.1932 

0.1932* 

^6 

0.25 

0.4 

0.35 

0.330 

d 

o 

0.328 

0.328* 

Q 





•H 



^rj 

1.0 

2.0 

1.50 

1.35 

-p 

05 CQ 
'HI 

1.31 

1.31 

^8 

0.01 

0.03 

0.015 

0.164 •' 

'H 

a -p 

•H 

0.0158 

0.0158 






-P fH 

- 19.20 


Xg 

15.0 

25.0 

20.0 

19.70 

PkO 

O 'H 

19.20 

Eximber 

of one' 

-d imensional 

9 


3 

2 

minimizat io ns 

for each 





unco ns t r a ined 

optimization 







tart , "• 


cO'/er rki'i lasua 
spar element 
















yo CHANG 



’/o CHANGE IN DESIGN VARIABLES 


THE CHANGE IN WEIGHT IS 
LESS THAN fA DUE TO 50% 
CHANGE IN X4 Xe-^g W; ' 







CHANGE IN DESIGN 
VARIABLES 


THE#eHANGE IN 8, IS LESS 
THAN 1%'Dui ■ TO- 50% CHANgI 
IN X, .Xcr .Xc 




■iill* 




li^ilKi|iPilH|S 




VARIATION OF MAXIMUM 






Vo CHANGE IN DESIGN VARIABLES' 


■ THE' CHANGE IN %.Root IS LESS 
‘ THAN. 1% DUE TO 50% CHANGE 
IN X 2 , 






illMlillil*®*** 


m: 











% CHANGE IN DESIGN 
VARIABLES 


THE "CHANGE INA^i IS LESS 
1% due TO 5G% ' CHANGE 


B.’S VA RI ATI 












2-.e rrecent apiroaoh is oapaile of solving the problem without 
:U> -indue :ti-dber of optinasation steps. M average of 21 


•ene 


i ,i::; oiiL'ic n 


a mimmization steps are required for optimzing a wing 
v.luii 9 decig a variables. 

(5) total computational time required for the optimization 

of one aircraft vang is about 45 minutes. Bor aircraft, the 
str-uctur:! v.ight saved can be converted directly into increased 
payload or indirectly into increased range. Hence the computatioml 
cost, ccmpaied to the reduction in weight obtained, can be seen to 
be very satrll. Horeover, designing a wing for strength, stability 
freo.uenoy (ani flutter) requiremente is certainly better than 
designing for strer^th alone and than satisfying the other requirements 
by mahi.^ necessary rnodmcations as is being done under «re present 

practice. 

thP behaviour of multiweb wing structures 
(6) In order to represent the behaviour 

the finite eluent idealisation using trian^ar membranes, reotan^ar 

Shear panels ari piifiointsd bars has been found to be simple and 

. . .1S.23 A reasonably aeourate analysis of the structure 

efficient , t 

only a fraction of the original 
can be performed even by oonsidermg ordy 

no +V, c+-«ir.tu-te Bae number of degrees of 
degrees of fJcaedom of the structure. 

^eedom in the static analysis can he reduced by ona-half by tsOcing 
the wi^ to be symmetric about itsrdddle plane and by using plate 
• aie order of the eigen value pioblem can 

fluxural assumptions. the o 

reduced by one-third by ^ploying the static condensatron 


further be 



technique 


The number of degrees of freedom in the flutter 


analysis can be reduced greatly by using the first few natural 
nodes as generalized coordinates instead of the nodal degrees 
of ff'eedom. 

(7 ) Sie total programme is mrritten in such a way that one can 
use ary of the requirements, either strength, stability, natural 
frequency, flutter and external shape constraints at a time or 
all of them put together. 

6.2 Recoimnendations 

( 1 ) Several improvements are possible at various stages of the 
present analysis. In the actual design, the number of degrees of 
freedom can be increased either by using backup storage or by 
making use of the theory of substructures. Instead of specifying 
a fixed-root condition for the wing the flexibility of the fuselage 
can also be considered in the analysis. 

(2) The possibilities of designing the wing under gust loads and 
landing loads can also be considered. As the aircraft has 
pass through various speed zones the unsteady aeix)dynamic loading 
can be refined by incorporating an aerodynamics package for subson c 
and transponic speed zones as well. Instead of designing the 
wing for a constant distributed load, the lift load resulting 
from a given flight condition can be used as design load. 
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(3) In order to compare the computational time a direct 
optimization method like Zoutendidk's feasible directions method 
can be aaployed for solving the optimization problaa. It is 

to be noted that this method does hot offer the flexibilily of 
using approximations in the analysis at several stages. 

(4) Since aerodynamic parameters like aspect mtio, sweepback 
angle are considered as design variables, it is a good caution 
to guard agaxnst drastic changes in the aerodynamic performance 

by placing constraints on the angle of attack, the lift and drag also 


(5 ) An obvious extension of the present work is to apply the 
automated minimum weight design procedure discussed in the present 
work to complete aircraft structures by treating the entire wing 
as a substructure and including fuselage, tails etc. ihis would 
require additional computer storage and different types of finite 
elements for proper representation of other structures of aircraft. 


* GOFOISJSIOHS - COJffHMJED 

8- It can also be concluded from the parametric studies 
that- in- tha initial .phase of minimum weight' design„of ^ s^pmetric, - 
: thin,- low a^ect ratio wing structures under bending inads^ the 
mass distributiDn of spars* and ribs can be ass^ 
thus reducing the number of variables -to , 

optimization. SMs will 

for computation. However, tbis requires mor^ -detailed study before 
arriving at a decision. 
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AFPEKDIX A 


Derivation of Element Stiftness and 
Mass Matrices 


The general procedure to be followed in deriving fjie 
stiffness and mass matrices of any finite element has been 
briefly outlined in Chapter 3. The detailed derivation of 
these ma-oi-ices from the assumed displacement or stress state will 
be press rrfced in this appendix for the three types of eloaents 
used in the idealization 


A- 1. Stiffness Matrix for a Triangular Membrane Element 

By assuming a linear variation of the displacement, the jnplane 
displacement functions u(x,y) and v(x,y) for the triangular plate 
shown in Eig, (A.i) can be presented as 



u(x,y) = d^x + d^y + d^ 

(a.i) 

and 

•v(x,y) = d^x + d^y + dg 

(A. 2) 


where the six aihitraiy coefficients d^, d^j.t.dg can be fourd 
from the Enplane displacements of the three vertices of the tri- 
angle. Using the boundary conditions 


and 


u(x^,y^) = U^ , v(x^,yp = 
uCxjry^ ) - = ■'^5 


(A.3) 
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in Bc[s, (iwl) and (a. 2 ) to evaluate the unknown coefficients, it 
can be shown that 



7/here 


= 2 A^ t (y3-y2) (x^- Xg) (y - y 2 ) 1 

= 2l^ Ky^-y^) (x- x^) -...(x^- x^) (y - y^)} (a. 5) 

= 2 ;^ — [ (y2“yT ) ) (y - yi ) ^ 


^12 ^14 ^16 ^21 “ ^25 " ^25 ° 

A^ 2 ^= Area of the triangle 123 = ^ ^^ 2 ^ 21 “ ^ 21 ^ 32 ^ 
and y^= x and y coordinates of node i in local 


coordinate system 



Bis. (a. 4) can now be used to find the relationship between the 

tot^ strains e , 1 and six -displacements U jV ,,.,7 , 

XX yy xy • ■ 

Differentiating these equations, we have 
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a?he stress-strain relation for the case of plane st3?oss is 


given by 


(a ^ 1 
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where E and v are Young's modulus and Poisson’ 
material. 


(A.7) 


s ratio of the 


The assumed displacement state satisfies the strain compa- 
tibility and the stress equilibrium equations within the element 
It can also be seen that the compatibility of displacements on 
two adjacent triangular elements with a common boundary is ensured. 
If JSqs, (a,4), (a, 6) and (A,?) are written in the matrix form 

u = [a] U 

t = [bl U , (A.8) 

and ? = Co] e 

the matrices [H , [b ] & to] can be identified from the 
equations (a,4), (a. 6) and (a. 7) respectively. 


(A. 6) 
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ilhe elanent stiffness matrix in local coordinate system can 
be obtained by evaluating the integral in (sj. {3.10)) , 

[k ] = I tb [ 0 ] [b] (a.9) 

Bie resulting matrix can be separated into two parts 

[k] - tkj + (bj 

where [k^] represents the stiffness matrix due to normal stresses, 

and tk^l represents the stiffness due to shearing stresses, Ohe 
tK> component matrices, as derived from m. (1.9) are given by 
Eqs. (a,ii) and (a.12) 


32 


^yiometric 


Et 

2 




^^123 

“'^^32^2 ^2 

2 




- ^32^31 ^^^2^31 

^31 



ft 1 « 

‘ tr 

'^^32^31 -"32^1 

"''=^31=^1 

2 

2 

^21: 


1 ^32^21 ■"'^^2^'^21 

j 

- =^3/21 



- -'^^ 32^21 ^2^21 

''^31='21 

“ ^1^21 

-vyj^si. 


:(A.ir) 
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Et 
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'("i'+v) 


^2 




”“^2^32 

2 

^32 



"^2^1 

^32^1 



^2^31 

”^32^31 •%#31 

2: 

%1 

2 

^21 

3^2^21 

“^32^21 “^1^21 

^31^21 

__ -^2^21 

^32^21 ^1^21 

”^3/21 

■ **21*^21' 


y. 


21 


(a, 12) 


vsbere t is the thic&nese the plate el^nenb. 
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Bl, (iUlO) gives the stiffness matrix in local eoordimte system 
(x,y). If the element is referred to some global coordinate 
system (x,y) as shown in Big. (A. 2), the new element stiffness 
matrix is given by 


[ h ] = [ X f [ k ] [ ^ ] 

9x9 9x6 6 x 6 6x9 

where [ X ] is a transformation matrix that relates the nodal 

-> 

displacements of the local coordinate system with those o± 
global system U by 


U = [ X ] 


6 

X 1 


6x9 

9 X 

1 

For the triangular membrane element 



“41 

41 

0 

0 


1^2 

“12 

m 

12 

0 

0 

[Xl = 

0 

0 

0 

1 

41 

^1 


0 

0 

0 

i 

12 

“^12 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


(A. 14) 


26 

the matrix [X] is given by 


0 

0 

0 

0 

0 

0 

0 

0 

“■41 

0 

0 

0 

“^12 

0 

0 

0 

0 


la, ^ 
41 

^^41 

0 

^12 

m „ 
12 

^2 


#iere 

*^41 " 
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“41 = % “i2S4'/'*43 

“41 * “l2'^14)/^43 ('‘•'6) 

^43 =t(^- + ( 73 - + G^- z^)^ - d^ 4 ] '''^^ 

A»2 Sftuivalent Mass Matrix for a Triangular Membrane Element 

Olae eq^iivalent matrix for the element in local coordinate 
system is given by (bi, (3.11 ))' 


[m 3 = / P [ a [a] dv 


(a. 17) 


where p is the density of the element and [ a ] has been defined 
in the Bg[, (ii,4). 

26 

It can be shown that for triangular plates undergoing 
essentially translational displacements (no bending involved), the 
mass matrix [m ] is invariant w-th respect to the coordinate 
transformation. Hence the element mass matrix for a traingular 
membrane element can be, obtained from Bq.s. (a,i 7) and (1.4) as 


C m ] =C m ] 
9x9 9 x9 


^^123^ 


0 0 1 0 0 1 0 0 
2 0 0 1 0 0 1 0 


0 0 2 0 
1 0 0 2 
0 1 0 0 


0 1 0 0 1 
0 0 1 0 0 


0 0 1 0 0 2 
1 0 0 1 0 0 


0 1 0 
0 0 1 
2 0 0 


0 1 0 0 1 0 0 2 
0 0 1 0 0 1 0 0 


(a. 18) 
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1*3 Stiffness Matrix for a Rectangular aiear Panel 


Itor the puipose of present analysis, liie trapezoidal shear 
panel is approximated by an equivalent rectangular dasar panel as 
shown in Pig, (a, 3). By assuming a linear stress distribution 
with in the rectangle 1234 


/ o \ + d-y 

I XX I 1 2 

Jvl = S-f 


(1.19) 


I ere d^,.,,d^ are constants, the displacement field u can be 


derived from the Hooke’s law 


0 = [c] e 


(A,20) 


Titiere [c] has been defined in the sq. (a. 7) 
Prom eqs. (a, 19) and (A,20), we obtain 


~ =~(d + dv- '’d-— V d .X ) 
3x E ^ 1 2^ 3 ^ 


=- 1(43 + 34 ,^- '> 4 ,- '’V) 


(A.21 ) 


(A. 22 ) 


^ 2(1+ v) (1.23) 

3y a X E 

solution of Eqs. (A.21) through (a. 23) can be given as 

2 d y^ 

u(x,y) « 4 td^x + d^xy - vd^x - vd^ 

2 2 

v(x.y) = i[-vV ^ 2 (i«)45X - 

dgX + dg ] (a, 25 ) 



99 


me unknown constants can be determined from 

the element displacements U and y. Ohe resxating strain 
displacement relation, specialized for -the case of a shear panel, 
can be expressed as 


e = [b] U 


(a.26) 



By using Bis, (a,?) and (A.27) in Bi. (a, 9) the stiffness matrix 
for a rectangular shear panel in local coordinate systmi can be 
evaluated and its final foarm is given by Bi. (A,28) 
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where G- is the shear modulus of the rectangular plate and ^ is 
the aspect ratio* 

3336 stiffness matrix of the element with respect to any 
global coordinate ^stem (xjy,z) can be obtained from 

rs]=[xf [k] [X] 

12 X12 12 X 8 8 x8 8x 12 

where the transfoimation matrix [x] is given by 



000000000 

000000000 

'23 “23 “23 ® ° ° ° ° ° 

^^2 “12 “12 ° ° ° ° ° ° (A.30) 

“ ° ° “23 “23 “23° ° ° 

° ° ° S2“l2“l2° ° ° 

0 0 0 0 0 0 Hjj Il>25 Ojj 

0 0 0 0 0 0 1^,2 “12 “12 



a.. 

ig 

= (x. - X. )/d. . 

1 3^' 13 



m . . 


(A.31) 


n. . 



ajcd 

d.. 

3.3 

= [(;.- ^ , f + (yj. yj^)^ A 
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It is to be Boted that the assumed stress distribution Bj, (a. 19) 
satisfies the stress equilibrium equations within the rectangle. 
Howerver, tiie resulting displacement distribution , Bqs, (a. 24) 
and (a, 25) violates the compatibility of boundary displacements 
on adjacent elements. 


Ipaped Mass Matrix for a Bectangular Shear Panel 


The equivalent mass matrix of a rectangular shear panel is 
given by the relation 


= / p [a ] ^ [ a ] dv (A.52) 

T 

where the matrix [a] can be obtained from B^s. (a. 24) and 
(a. 25). However, the lumped mass matrix of the equivalent 
rectangular plate shown in Pig- (A. 5) is used for the trapezoidal 
shear panel also as an approxiiaation. It is given by 




whe 3 ?e 


a is the depth of the element and b is the span* 


A.5 Stiffness Matrix for a Pin-Jointed Bar Element. 

A pior- jointed bar is a on^ dimensional element for vfcich the 
assumed displacement can be taken as [see Pig. (a.4) ] 

U. 


u(x) = [a] U = [(i-j) -IjJ ^ 




(a* 34) 


2.J 


where ^ is the of the bar. 

9u 


and 



(A. 35) 




(A.36) 


The element stiffness matrix in local coordinate system is given by 

“ 1 


[k] = / [b] ® [c] [b] dv = 


x2 


P { 


-1 


1 


(A* 37) 


where A is ihe cross-sectional area of the bar# 

Here, the nodal displacements in local and global coordinate 
^sterns are related by 

ft 

1 




U 


[X] U 5 


^ m n 0 0 0 

12 “l2 12 

0 0 0 ^-,2 ^-, 2^2 




tJ. 








(A. 38) 
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•where given by E^s. (a.3i) 

She stiffness matrix of the bar with respect to the global 
coordinate systaa is given by 


[k] 

6 x6 


= [X] 

6 x2 


T 


[k] [ X] 

2x2 2x6 


(A.39) 


A, 6 a3,uiva3.ent Mass Matrix for a pin. Jointed Bar ELement 

Prom Fig. (a, 4 ), the displacements at ai^r point along the length 
of the bar can be obtained as 


u(x)'| 


0 

0 

x/^ 

0 

0 

•| v(x) j = 

0 

1-x/®' 0 

0 

x/l 

0 

^ w(x) j 

0 

0 

1-x/ & 

0 

0 

x/ Z 




V. 




u. 






(a. 40) 


or 


u = [a] U 

3x'| 3x6 6x1 


(A.41) 


fhas the matrix [a ] can be identified as 

Wa 0 0 Va 0 0 

0 Wii 0 0 -x/i 0 

0 0 0 0 'x/a 


Ea(x)3= 


(a. 42) 
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Assuming the cross-sectional area of the bar A to be constant, 
the equivalent mass matrix is given by 

2 0 


[ m] = [m ] = /p [a]^{a] dv = 
6x6 6x6 - 


0 

0 

1 

0 

L? 


2 

0 

0 

1 

0 


0 

0 

2 
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0 

1 


1 

0 

0 
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0 

c 


0 

1 
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0 

2 

0 
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0 

1 

0 

0 

2 


tA.43) 
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APPEiaiX B 


asses ipnoN op mr rn-TriA-:':’ 


!I?he ccHttputer piogramme is written in it»rtriin iv lor 

TTOff 7044 computer. ®ie progranine consists of irt subitjutinew wKlc'r. 
are called by MAIN during optimization sti/ige. Tiiia pro/rramc c'ln 
be used to analyze or design argr win^j structure th;tt could be 
represented by cover skins, v/ebs aM stringers* 'Pue 
calculates deflections, stresses, mturtal i’requcncies, flutter Jipeeii® 
and flutter frequencies. In design mode the prOf-rrrjanie can be used 
to find the minhmM weight of wing structure satisfying the otrengthf 
stability, frequency, flutter and geometric comtitiinta* 

The input data consists of the following informatioiM 

( 1 ) Details of the finite element modeling and node numbering acherae. 

(2) Material psroperties of the wing 

(3) Details of the payload and flight conditions 

(4 ) A feasible starting design vector 

(5 ) Convergence criteria 

'Bie following output can be obtained: 

( 1 ) .'HI the input information 

(2) An analysis of the wing, if required 

(3) She intenaediate azd final ^resulte of anrilysis and optiaiKi» 
tlon in the design mode* 
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B« 1 Purpose of the sabrou tines 


.( 1 ) 


It xii 5 )l€ments Davidon-Pletcher-Powell variable metric 
method of unconstrained minimization. 


(2) PUl'I 


(3) HJII1 


(4) 

OBJ 

( 5 ) 

XTCORD 

(6) 

ADATA 

( 7 ) 

STORE 


(8) Slip 


(9) PUGE 


It evaluates the penally function and its gradient, if 
the change in design vector is small, using linear 
appioximation, otherwise exact analysis. 

Constraints are formulated. Design vector is checked 
whether it is feasible. The penalty function value is 
evaluated. 

It evaluates the weight of wing 

It evaluates the dimensions of wing and correlates the 
design vector to the gauge parameters of the elements 
It i^ads all the input information 

It assembles stiffness matilx [SC], mass matrix [31] , 
aerodynamic matrices [SA] , [SB] , [SC] and [SD] , Then 
it performs static analysis, reduces the stiffress and 
mass matrices using static condensation and obtains the 
natural frequencies and the associated vectors. Afterwords it 
calls PDUTER to perform the double eigenvalue analysis 
It determines the type of the element and calls the 
appropriate subroutine to get elemental stiffness and 
mass matrices 

It gives elemental stiffness and mass matrices for the 
pin— jointed bar element. Also using the displacement 
information, it calculates the stresses in the specified 


elements. 
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(1 0 ) SHEAR It performs the same functions as iPLlirGE’ for triangular 

membrane elements, besides giving the elemental aero- 
dynamic matrices 

(11) IiCfWGI It performs the same functions as ’EIiH'GE’ for rectangular 

shear panels 

(12) AGROSB It multiplies two real matrices 

(15) ASRIfSB It also multiplies two real matrices, however it give^f 

[A®]x[E] 

(14) JOEDiJff It finds the inverse of a real matrix 
( 15 ) MIQ It solves the eigenvalue problem using power method 

(16) SOLVE It solves [K] ? =P 

(17) EIUIEK It Impi ements the double eigenvalue analysis of flutter 

(18) GCMABI It finds the value of a complex determinant 

Some of the main features of these subroutines axe presented in 
the form of flow diagram in the next section. The flow charts for 
other subroutines which use standard algorithms are not included in 
this flow diagram. 



B,2 Slow Diagram 


( T ) MAIS 
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(2) DJ!P 
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(3) HJIT 
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(4) iUJT 1 







(5 ) S20EE 



Obtain dimensions of wing 
from XYCORf 
Initialize [SK] 


Il^ 

64 

DO KOUDT-I 4 : 




Assemble [SCl 

and, store it as [affi] 

obtain displacements 
from SOITE 


IS 

KOIJITT = 1 ? 


^ IS ^ 
KOURT = 2 ? 


Perform stress ai^ysis 
Restore [SS] & R1 
Initialize [3M] 
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Continue 


IS 

KOURT = 3 


Assemble [3)1] 


Assemble [SA] , [SB], [SC], [SD] 


Perform static condensation 
obtain eigen solution from lEBQ 


Obtain Flutter solution 
from FIUfER 















91881 


A£ ,r ~ K i ^ M 



